package it.uniroma2.dtk.op.convolution;

import java.util.Arrays;

import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;

import it.uniroma2.dtk.op.TransformAndCompose;
import it.uniroma2.util.math.Complex;
import it.uniroma2.util.math.FFT;
import it.uniroma2.util.vector.VectorProvider;

public class CircularConvolution extends TransformAndCompose {

	protected FFT fft = null;		// Object used to optionally compute the Fast Fourier Transform for the circular convolutions
	protected int spaceGrade = 0;	// log_2 of the vector dimension, where an integer, for using the FFT
	
	protected DoubleFFT_1D fft2 = null;
	
	@Override
	public double[] op(double[] x, double[] y) {
		if (x == null)
			return y;
		else if (y == null)
			return x;

		if (fft2 != null) return circularConvolutionFFT2(x, y);
		else if (fft != null) return circularConvolutionFFT(x, y);
		else return circularConvolutionBasic(x, y);
	}
	
	protected double[] circularConvolutionBasic(double[] firstVector, double[] secondVector) {
		int size = firstVector.length;
		double[] result = new double[size];
		Arrays.fill(result, 0);
		for (int i=0; i<size; i++)
			for (int j=0; j<size; j++)
				result[i] += firstVector[j] * secondVector[(i-j)<0 ? (i-j+size) : i-j];
		return result;
	}
	
	protected double[] circularConvolutionFFT(double[] firstVector, double[] secondVector) {
	return Complex.extract(fft.real_cconvolve(
					Complex.generate(firstVector),
					Complex.generate(secondVector),
					spaceGrade));
	}
	
	
	protected double[] circularConvolutionFFT2(double[] a, double[] b) {
		int N = a.length;
		double [] first = new double[2*N];
		for (int i=0; i < N ; i++) first[i] = a[i];
		double [] second = new double[2*N];
		for (int i=0; i < N ; i++) second[i] = b[i];
		fft2.realForwardFull(first);
		fft2.realForwardFull(second);
	    double [] c = new double[2*N];
	    for (int i = 0; i < N; i++) {
	        c[2*i] = first[2*i]*second[2*i] - first[2*i+1]*second[2*i+1];
	        c[2*i+1] = first[2*i]*second[2*i+1] + first[2*i+1]*second[2*i];
	    }
	    // compute inverse FFT
	    fft2.complexInverse(c, true);
	    double [] out = new double[N];
	    for (int i=0;i<N;i++) out[i] = c[2*i];
	    return out;
	}
	
	
	
	@Override
	public void initialize(VectorProvider vp) throws Exception {
		super.initialize(vp);
		int exp = log2(vp.getVectorSize());
		if (exp < 0) System.out.println("Using the slow circular convolution as vector size is not power of 2");
		else {
			if (System.getProperties().containsKey("fft1")) {
				System.out.println("Fast CircularConvolution in use");
				fft = new FFT();
				fft.initialize(exp);
			} else {
				System.out.println("Fast CircularConvolution (with threads) in use");
				fft2 = new DoubleFFT_1D(vp.getVectorSize());
			}
			spaceGrade = exp;
		}
	}
	
	protected int log2(int arg) {
		int exp = 0;
		while(true) {
			double pow = Math.pow(2, exp); 
			if (pow == arg)
				return exp;
			else if (pow > arg)
				return -1;
			exp++;
		}
	}

}
